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COMPUTING PLANAR AND SPHERICAL CHOREOGRAPHIES 

HADRIEN MONTANELLI* AND NIKOLA I. GUSHTEROVt 


Abstract. An algorithm is presented for numerical computation of choreographies in the plane in a Newtonian 
potential and on the sphere in a cotangent potential. It is based on stereographic projection, approximation by 
trigonometric polynomials, and quasi-Newton and Newton optimization methods with exact gradient and exact 
Hessian matrix. New choreographies on the sphere are presented. 
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1. Introduction. Choreographies are periodic solutions of the n-body problem, n > 2, in 
which the bodies share a common orbit and are uniformly spread along it. The study of chore¬ 
ographies with unit masses in the plane in a Newtonian potential is an old one. For the two-body 
problem, the only choreography is a circle, and for the three-body problem, the first one was found 
by Lagrange in 1772 m , and is also a circle. The second choreography of the three-body problem 
with unit masses, a figure-eight, was discovered numerically more than two centuries later by Moore 
in 1993 [15) . while Chenciner and Montgomery proved the existence of this class of orbits a few 
years later [4]. In the early 2000s, many new choreographies for various n and with unit masses 
were found by Simo using a combination of different numerical methods |21| . In 2008, Chen proved 
the existence of infinitely many choreographies of the three-body problem with a certain topological 
type for various choices of masses 0. 

The situation is different on the sphere 0 . Although there is a growing interest in the n-body 
problem on the sphere (and in other spaces of constant curvature) in a cotangent potential (T, [2j 
0 Hinnsi m ini, the only non-circular choreographies with unit masses found so far are for the 
two-body problem pB 

We present in this paper an algorithm to compute planar and spherical choreographies with 
unit masses to high accuracy. We have found many new choreographies on the sphere of radius R 
in a cotangent potential for various n > 2. These are curved versions of the planar choreographies 
found by Simo and, as R -A oc, they converge to the planar ones at a rate proportional to the 
curvature 1/R 2 . 

2. Planar choreographies. Let Zj(t ) G C, 0 < j < n — 1, denote the positions of n bodies 
with unit mass in the complex plane. The planar n-body problem describes the motion of these 
bodies under the action of Newton’s law of gravitation, through the nonlinear coupled system of 
ODEs 


n —1 


m - e 


2 = 0 


z i(t) ~ z j(t) 

z i(t ) - zj(t) | 3 


= 0 , 


0 < j < n — 1. 


( 2 . 1 ) 
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1 Throughout this paper, the term “sphere” will always refer to the 2-sphere. 

2 Diacu proved the existence of choreographies with unit masses for the two- and three-body problems on the 
3-sphere in [6]. 
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We are interested in periodic solutions of m in which the bodies share a single orbit and are 
uniformly spread along it, that is, solutions Zj(t ) such that 


Zj(t) = q(t + o < j < n — 1, 


( 2 . 2 ) 


for some 27r-periodic function q : [0, 27r] -A C. Such solutions were named choreographies by Simo, 
the n bodies being “seen to dance in a somewhat complicated way” m The period can be chosen 
equal to 27r because if q{t) is a T-periodic solution of (|2.ip . then A _2 / 3 g(At), A = T/( 27t), is a 
27r-periodic one. It has been well known since Poincare [18, 19] that the principle of least action , 
first introduced by Maupertuis in 1744 m, can be used to characterize periodic solutions of (£D): 
choreographies (IQi are minima of the action functional , or simply action , defined as the integral 
over one period of the kinetic minus the potential energy, 


with kinetic energy 


A = 


p2tt 

/ (K(t) - U(t)) dt, 

Jo 


1 n —1 1 n— 1 0 . 


i=0 


7=0 


and potential energy 


(2.3) 


(2.4) 


n-lj-l 

u*) = -EEh(*uu 

j=02=0 


1-1 


n—17—1 


ft —-l j ■*- o • 




J=02=0 


2frj \ 
n / 


(2.5) 


Note that the action (|2.3p depends on q(t) via U(t) and on q'(t) via K(t). Since the integral of 
(|2.4p does not depend on j and the integral of ([2.5D only depends on i — j, the action functional 
can be rewritten 



( 2 . 6 ) 


Planar choreographies correspond to functions q(t) which minimize ([2.6D . 

We are also interested in solutions of (E3D in which the bodies share a single orbit q(t) that is 
rotating with angular velocity uj relative to an inertial reference frame, i.e., 

Zj(t) = e luJt q(t + , 0 < j < n — 1. (2.7) 

Choreographies of the form (1277)1 are said to be relative , as opposed to the absolute choreographies 
(12.21) . The action associated with relative planar choreographies is 



+ iuoq(t)\ 2 dt - 


n 

2 


**-— 1 2 tt 




) 


-1 

dt. 


( 2 . 8 ) 


Note that (|2.6p is the special case of (|2.8p with uj = 0. 
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3. Computing planar choreographies. Our method for computing planar choreographies 
is based on the minimization of the action (|2.8|) and uses two key ingredients: 

Ingredient 1. Trigonometric interpolation. The function q(t) is represented by its trigonometric 
interpolant in the exp {ikt) basis. The optimization variables are the real and imaginary parts of 
its Fourier coefficients. The action is computed with the exponentially accurate trapezoidal rule. 

Ingredient 2. Closed-form expressions for the gradient and the Hessian. Formulas for the 
gradient and the Hessian matrix of the action (|2.8|) with respect to the optimization variables are 
derived explicitly and used in the optimization algorithms. 

The numerical optimization of the action is in two steps: 

Step 1. Quasi-Newton optimization methods. Numerical optimization methods with the exact 
gradient and based on approximations of the Hessian are employed with a small number of opti¬ 
mization variables. The accuracy of the solution at this stage is from one to five digits. This step 
is computationally very cheap. 

Step 2. Newton’s method. Once an approximation to a choreography has been computed via 
a quasi-Newton method, one can improve the accuracy to typically ten digits with a few steps of 
Newton’s method with exact Hessian, and a larger number of optimization variables. This step is 
computationally more expensive. 

Let us start with a few words about the first ingredient. The approach used by Simo m 
is to decompose the function q(t) into real and imaginary parts, and to represent each of them 
by a trigonometric interpolant in the sin (kt) and cos(kt) basis. In this paper, we use instead a 
trigonometric interpolant of the function q(t) itself in the exp (ikt) basis. For an odd number AT, let 
{tj = 27rj/N}, 0 < j < TV —1, denote N equispaced points in [0, 27r) and {qj = q(tj )}, 0 < j < AT —1, 
the (complex) values of q(t) at the tj’ s. The trigonometric interpolant pN(t) of q(t) at these points 
is defined by 


with Fourier coefficients 


PN{t) = T, c kC tkt , t 6 [0, 2n], 

k - N-l 


N—l 




AT- 1 


3=0 


(3.1) 


(3.2) 


The trigonometric interpolant problem goes back at least to the young Gauss’s calculations of the 
orbit of the asteroid Ceres in 1801—it seems that planetary orbits and trigonometric interpolation 
share a long and on-going relationship. Throughout this paper, the number of grid points AT 
will always be odd. All our results have analogues for AT even, but the formulas are different, 
and little would be gained by writing everything twice. If we replace q(t) by its trigonometric 
interpolant <pni) - (poi) with Cfc = u k + ivk, the action (|2.8|) becomes a function of the 2AT real 
variables {uk,Vk}, \k\ < (AT — l)/2. We are looking for solutions q(t) without collisions. The 
integrands in ([2.8D are therefore analytic and the trapezoidal rule converges exponentially [22]. We 
use Chebfun v5.2. im to compute trigonometric interpolants. Chebfun is an open-source package, 
MATLAB-based, for computing with functions to 16-digit accuracy. Its recent extension to periodic 
functions [©J provides a very convenient framework for working with closed curves in the complex 
plane. 
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Let us now say more about the second ingredient. The exact gradient and exact Hessian are 
derived in Appendix A. The gradient can be computed in 0{nN 2 ) operations, while the computation 
of the Hessian requires 0(nN 3 ) operations. 

The numerical optimization of the action, in MATLAB R2015b, is carried out in two steps. 
First, we apply a quasi-Newton method m Chapter 6] using the exact gradient, and with a small 
number of Fourier coefficients, N = 55 and 75 in our experiments. Quasi-Newton methods are 
based on the approximation of the Hessian matrix (or its inverse) using rank-one or -two updates 
specified by gradient evaluations. In MATLAB, the fminunc command implements various quasi- 
Newton methods and, among them, we choose the BFGS algorithm m ■ We take O(N) iterations 
of the BFGS algorithm. The cost of this first step is thus 0(nN 3 ) since, at each iteration, BFGS 
computes the gradient in 0(nN 2 ) operations and matrix-vector products in 0(N 2 ) operations. 
Second, we perform a small number s of iterations of an approximate Newton method with exact 
Hessian and M Fourier coefficients, M > N . The starting point of the approximate Newton 
method is the output of the BFGS algorithm, padded with M — N zeros. An exact Newton 
method would have a 0(snM 3 ) cost, since, at each iteration, it requires the computation of the 
exact Hessian (0(nM 3 ) operations) and the solution a linear system (0(M 3 ) operations). To 
reduce this cost, we use an approximate Newton method. We compute the exact Hessian at the 
first iteration only, and compute its LDL T decomposition (0(M 3 ) operations, a generalization of 
Cholesky decomposition for symmetric matrices which are not positive definitive, typically half 
as expensive as LU factorization). The first iteration has thus a 0(nM 3 ) cost. The subsequent 
iterations of Newton’s method do not recompute the Hessian but use this factorization instead. The 
solution of the linear system can then be computed in 0(M 2 ) operations at each iteration. The total 
cost of this second step is then 0(nM 3 ), and the total cost of the optimization is 0(n(N 3 + M 3 )). 

Let us add four comments about this optimization process. First, since the initial guess of 
Newton’s method—the output of BFGS—is a good approximation of a choreography, the Hessian 
matrix does not vary significantly from one iteration to another. As a consequence, using the LDL T 
factorization of the Hessian of the first iterate at each iteration does not affect the convergence very 
much. Second, at a minimum of the action, i.e., a choreography, the Hessian is positive definite, so 
we could in principle use the Cholesky decomposition instead of the LDL T decomposition. However, 
in practice, because the Hessian is computed at an approximation of a choreography, it often has 
some small negative eigenvalues. Third, we do not use Newton’s method with exact Hessian from 
the beginning because it only converges for initial guesses close enough to the solution. Fourth, 
another option would be to use M coefficients with the quasi-Newton method directly. However, 
we found in practice that BFGS with M coefficients typically achieves an accuracy of 6 digits at 
most, while Newton’s method achieves an accuracy of 10 digits. 

For both steps of the optimization, the accuracy is defined as the the 2-norm of the residual of 
(|2.1|) divided by the 2-norm of the solution (relative 2-norm). The residual is computed in Chebfun 
with the chebop class m, the Chebfun automatic solver of differential equations. We also check 
that the 2-norm of the gradient divided by the 2-norm of the gradient of the initial guess (relative 
2-norm) is close to zero, and that the Fourier coefficients of the solution decay to sufficiently small 
values. 

The famous figure-eight, with action A « 24.371926 |4], is plotted in Figure l3Tl It is obtained 
by running the code of Figure lT3l The code uses the actiongradeval and gradhesseval functions, 
which compute the action and the gradient, and the gradient and the Hessian; the codes are avail¬ 
able online at the first author’s GitHub web-page (http://github.com/Hadrien-Montanelli). 
Table 13.11 shows some numbers pertaining to the computation of the figure-eight, including the 
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Fig. 3.1. The famous figure-eight solution of the three-body problem, obtained by running the code of Figure 13.31 
The dots show the bodies at time t — 0. The action of the resulting choreography, A = 24.371926476242812, agrees 
with the 8 digits given in 0- 


Fourier coefficients 
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O BFGS + two iterations of Newton 
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Fig. 3.2. Absolute values of the Fourier coefficients of the figure-eight of Figure [3~Tt obtained by BFGS (red 
dots) and BFGS followed by two iterations of Newton’s method (black circles). 


relative 2-norms of the solution and the gradient, and the amplitude of the smallest (numerically 
nonzero) Fourier coefficient. After 51 iterations of the BFGS algorithm, the solution is accurate to 
five digits, and two iterations of Newton’s method gives six extra correct digits. The solution found 
by BFGS and BFGS plus Newton look the same to the eye; if they were plotted on the same graph, 
they would be perfectly superimposed. The difference is visible in coefficient space. We plot the 
Fourier coefficients of the two solutions in Figure [3721 Since choreographies are analytic functions, 
the Fourier coefficients decay geometrically [23] Theorem 4.1]. The solution obtained by the BFGS 
algorithm only uses 55 coefficients, that is, wavenumbers \k\ < 27, and the Fourier coefficients 
decay to about 10 -8 . The solution obtained by Newton’s method uses 145 Fourier coefficients, i.e., 
\k\ < 72, which decay to machine precision. 

All the planar absolute choreographies of the five-body problem found by Simo El Figure 2] 
can be computed with this algorithm. We plot six of them in Figure l3~4l Tables lT2l and EPl show 
some numbers pertaining to their computation. The BFGS algorithm leads to results accurate to 
a few digits, and two to six iterations of Newton’s method lead to about ten digits of accuracy. 
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°/ 0 Initial guess: 

n = 3; N = 55; M = 145; 

qO = chebfun(@(t)cos(t)+li*sin(2*t), [0 2*pi],N,'trig'); 
cO = trigcoeffs(qO); 

°/ 0 BFGS algorithm: 

options = optimoptions('fminunc'); 
options.GradObj = 'on'; 
options.Algorithm = 'quasi-newton'; 
options.HessUpdate = 'bfgs'; 

c = fminunc(@(x)actiongradeval(x,n),[real(cO);imag(cO)],options); 

°/ 0 Two iterations of an approximate Newton method: 

c = [zeros((M-N)/2,1);c(1:N);zeros(M-N,1);c(N+l:end);zeros((M-N)/2,1)]; 
mid = 1 + floor(M/2); 

[G, H] = gradhesseval(c,n); [L, D] = ldl(H); 
for k = 1:2 

s = L'\(D\(L\(-G))); 

cnew = c + [s(1:mid-l);0;s(mid:M+mid-2);0;s(M+mid-l:end)]; c = cnew; 

G = gradhesseval(c,n); 
end 

°/ 0 Reconstruct solution: 

c = c(1:M) + li*c(M+l:2*M); 
q = chebfun(c, [0 2*pi],'coeffs','trig'); 

Fig. 3.3. MATLAB code to compute the figure-eight of Fiaure \3A\ This code gives correct results to 11 digits 
of accuracy in less than half a second on a 2.7 GHz Intel i7 machine. 



-1 0 1 


-1 0 1 


-1 0 1 





Fig. 3.4. Some absolute planar choreographies of the five-body problem. These correspond to choreographies 1 
(top-left comer), 2 (top center), 4 (top-right corner), 7 (bottom-left corner), 16 (bottom center) and 18 (bottom-right 
corner) of H3 Figure 2]. The dots show the bodies at t = 0. 
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BFGS 

Newton 

Action 

24.371926476245442 

24.371926476242809 

Number of coefficients 

55 

145 

Computer time (s) 

0.19 

0.27 

Number of iterations 

51 

2 

Relative 2-norm of the gradient 

2.86e-07 

3.90e-16 

Smallest coefficient 

3.17e-08 

2.09e-18 

Relative 2-norm of the residual 

2.06e-05 

2.24e-ll 


Table 3.1 

Computation of the figure-eight choreography of Figure 13.11 



1 

2 

4 

7 

16 

18 

Action 

68.8516 

71.3312 

77.1588 

88.4397 

109.6361 

119.3168 

Number of coefficients 

75 

75 

75 

75 

75 

75 

Computer time (s) 

0.89 

0.63 

0.59 

0.62 

1.73 

0.83 

Number of iterations 

98 

68 

69 

136 

373 

166 

Relative 2-norm of the gradient 

3.27e-08 

6.65e-08 

3.36e-07 

4.85e-09 

3.21e-10 

2.89e-09 

Smallest coefficient 

1.28e-05 

6.13e-09 

5.89e-06 

8.12e-06 

3.61e-05 

6.52e-05 

Relative 2-norm of the residual 

3.51e-02 

1.07e-05 

Table 3.2 

6.65e-03 

2.23e-02 

2.75e-01 

3.55e-01 

Computation of the absolute 

: planar choreographies of Figure 13.41 

with the BFGS algorithm. 




1 

2 

4 

7 

16 

18 

Action 

68.8516 

71.3312 

77.1588 

88.4397 

109.6361 

119.3184 

Number of coefficients 

335 

145 

245 

285 

445 

455 

Computer time (s) 

3.49 

0.48 

1.68 

2.33 

7.26 

7.85 

Number of iterations 

4 

2 

3 

4 

5 

6 

Relative 2-norm of the gradient 

6.93e-16 

1.90e-16 

1.71e-13 

1.53e-16 

3.74e-14 

3.79e-14 

Smallest coefficient 

3.43e-15 

8.81e-16 

8.00e-17 

2.05e-14 

1.26e-16 

1.43e-14 

Relative 2-norm of the residual 

8.29e-10 

5.66e-ll 

7.27e-10 

8.31e-10 

1.51e-09 

2.62e-09 


Table 3.3 

Computation of the absolute planar choreographies of Figure E3] with Newton’s method using the outputs of 
BFGS as initial guesses. 


The same method can be used to compute relative choreographies to high accuracy. We plot 
three relative planar choreographies of the seven-body problem in Figure [3751 

An interactive tool to compute choreographies with MATLAB and Chebfun is available at the 
web-page previously given. The code, choreo, finds choreographies starting with hand-drawn initial 
guesses. It is easy to use, fast and enjoyable—the reader is highly encouraged to try it! 

Let us conclude this section with a few words about the number of choreographies for a given 
n. This number is not known, but there is an interesting result, due to Simo m Proposition 5.1], 
about the (smaller) number of choreographies that consist of a concatenation of “bubbles,” such as 
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Fig. 3.5. Relative planar choreographies of the seven-body problem with angular velocity 2.8 (left), —2.9 (center) 
and 2.31 (right). The dots show the bodies at time t = 0. They can be computed to ten digits of accuracy with about 
four hundred Fourier coefficients. 


choreographies 1, 2 and 4 of Figure l3~4l For n > 3, there are 2 n 3 + 2^ n 3 )/ 2 J such choreographies. 

4. Spherical choreographies. Let Xj(t) G M 3 , 0 < j < n — 1, denote the Cartesian coor¬ 
dinates of n bodies with unit mass on the sphere = {X G M 3 , ||X|| = i?}, where || • || is the 
Euclidean norm in M 3 . The n-body problem on the sphere in a cotangent potential describes the 
motion of these bodies via the n coupled nonlinear ODEs 


' ^ [R* - (Xi(t) ■ Xj^f 2 


+ R- 2 \\X' j {t)\\ X j {t) = 0, 0 < j < n — 1. (4.1) 


i =0 
A 7 


See [Hj for details about the derivation of these equations. Note that the potential associated with 
m is no longer the Newtonian potential (12311 . It is a cotangent potential, a generalization of 
the Newtonian potential on the sphere, and dates back to the 1820’s with the work of Bolyai and 
Lobachevsky. The reader can find a detailed history of the problem in Diacu’s 2012 book on relative 
equilibria [5]. 

We are looking for periodic solutions of 63D moving along the same orbit, i.e., solutions Xj (t) 
such that 

Xj(t)=Q(t+^, 0 < j < n — 1, (4.2) 

for some 27r-periodic function Q : [0,27r] -a C M 3 . Again, the period can be chosen equal to 
27r because if Q(t) is a T-periodic solution of (|4.1j) on the sphere of radius i?, then A _2 / 3 Q(At), 
A = T/(27 t), is a 27r-periodic one on the sphere of radius A ~ 2 ^R. We call these solutions spherical 
choreographies. They are minima of the action associated with (|4.1[h defined again as the integral 
over one period of the kinetic minus the potential energy, with kinetic energy 


and potential energy 


1 n —1 1 n —1 

m = - 2 Z\\mt = 2T,M*+ 


j =0 


3=0 


27TJf\ 
n J 


n— 1 j —1 


f/ w = -^EE cot 




3=0 i =0 


(4.3) 


R 


(4.4) 
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where 


d(Xi(t),Xj(t)) = R arccos ^^ — (4.5) 

is the great-circle distance between Xi(t) and Xj(t ) on The potential ([4.41) is the cotangent of 
the (rescaled) distance on the sphere. Using the trigonometric identity cot(arccos(x)) = x/y/1 — x 2 , 
the potential energy can be rewritten 


n— 1j—1 


-| -L .7 


*<(*)•*,•(*) 




(4.6) 


The action is then given by 



Q(t) ~ Qjt + 

V /^-(Q(t).Q(t+2£l)) 2 


dt. 


(4.7) 


Spherical choreographies correspond to functions Q(t) which minimize (|4.7|) . Note that since the 
cotangent potential (|4.4|) is singular not only when the distance between two bodies is zero but also 
for antipodal configurations, we are looking for solutions that stay in a single hemisphere. See m 
for more details about the singularities of the n-body problem in a cotangent potential. 

As in the plane, we are also interested in solutions of (|4.1|) in which the bodies share a single 
orbit Q(t) that is rotating with angular velocity uj along the 2 )-axis relative to an inertial reference 
frame, i.e., 


Xj(t) = 


cos (cut) 
sin(co’t) 
0 


— sin (cut) 
cos (cut) 
0 



0 < j < n — 1. 


(4.8) 


Let R u (t) denote the rotation matrix in (14.8)1 . The action associated with relative spherical chore¬ 
ographies is 


A = 


n f 2n 

2 Jo 


\R u {t)Q , {t) + R! u (t)Q(t)\\ 2 dt- 


n—1 2tt 

*5i 


Q(t) ■ Q(t + 


R*-(Q(t)-Q{t+ 2 -£)Y 


: dt. (4.9) 


As in the plane, (in is the special case of (|4.9|) with uj = 0. 

5. Computing spherical choreographies. Our method for computing spherical choreogra¬ 
phies is based on stereographic projection and on the algorithm described in Section 3. Points 
X = (xi,X2,£3) t on the sphere are mapped to points z = Pr(X ) in the plane C via 


* = Pr(X) 


Rx i + iRx 2 
R- x 3 


(5.1) 


_r 2 TjT |z|2 (2-R 2 Re(z), 2R 2 lm(z), - R 3 + i?|*| 2 ) T 


The inverse mapping is given by 
X = P^(z) = 


(5.2) 














10 


COMPUTING PLANAR AND SPHERICAL CHOREOGRAPHIES 


The Euclidean distance d(X,Y) = \\X — Y\\ between two points on the sphere is transformed into 
the distance d(z,£) between their projections z = Pr(X ) and £ = Pr(Y ) defined by 

d(z, 0 = 2R 2 \z~Z\ (5.3) 

vWUW+f) 

and the great-circle distance (1T51) into 

d(z,£) = 2R arcsin (5.4) 

The complex plane endowed with the distance 63} is called the spherical plane. Let q{t) = PR^Qit)) 
denote the projection of the curve Q(t) onto C, and 


Zj(t) = P R (X j (t)) = P R {Q(t+^j) =q(t + ^, 0 <j <n—l, 

the projections of the n bodies Xj(t). The action (|4.9|) can be then reformulated as 
4 = n f 2 ” f 2R?\q'(t)+iu>q(t)\ \ 2 ^ , f 2 ” 2 & - Ptf) 2 

2 Jo \ R 2 + \q(t) | 2 7 2R ^Jo Dj{t)y/AR 2 -Dj{t ) 2 


(5.5) 


(5.6) 


with Dj(t) = d(q(t). q(t + AT ■ Perez-Chavela and Reyes-Victoria |1T] Theorem 2.3] showed the 
equivalence of the formulations (|4.9D and (|5.6D with cj = 0. 

Once the problem is reformulated in the spherical plane, we apply the two key ingredients 
and the two steps described in Section 3. The function q(t) is approximated by its trigonometric 
interpolant (|3.ip ~ (|3.2|) at N points, the action (|5.6| ) becomes a function of the real and imaginary 
parts of the Fourier coefficients, and is computed with the exponentially accurate trapezoidal rule. 
Formulas for the gradient and the Hessian matrix of the action (|5.6 | ) are derived in Appendix B 
and used in the optimization algorithms. Again, the computation of the gradient costs 0(nN 2 ) 
while the computation of the Hessian requires 0(nN 3 ) operations. For the optimization, we use the 
same strategy: BFGS algorithm with exact gradient and a small number of variables, followed by a 
few steps of an approximate Newton method with exact Hessian and a larger number of variables. 
As in the plane, at convergence, we check that the norm of the gradient of the action is close to 
zero, the Fourier coefficients of the solution decay to sufficiently small values, and the solution 
satisfies equation 63} projected into the plane. The latter was first given by Perez-Chavela and 
Reyes-Victoria in 2012 on Lemma 2.1], and can be written as 


m 


R 2 + \zM 2 


4 R Pj t i(t) 


0 < j < n — 1, 


(5.7) 


where A j(t) = 4 R A / ( R 2 + \zj(t)\ 2 ) 2 is the conformal factor that appears in the kinetic part of (|5.6|h 
while Pj,i(t) and Oj,i(t) are defined by 


Pj,i(t) = [R 2 + \ Zj (t )| 2 ] [. R 2 + \ Zi (t)\ 2 ] 2 [R 2 + Zi(t)zj{t)\ [. Zi (t ) - Zj(t)], (5.8) 

%,T) = - [2 R 2 z j {t)z i {t) + 2R 2 Zi(t)zj(t) + (| zj(t )| 2 - R 2 )(\ Zi (t)\ 2 ~ R 2 T q] 

+ [R 2 + \ Z j(t)\ 2 ] 2 [R 2 + |zj(i)| 2 ] 2 . 
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Fig. 5.1. A circular choreography of the spherical three-body problem on the sphere of radius R = 1 (left) and 
its projection in the plane (right). All the circles of radius r, 0 < r < R, are spherical choreographies for any R ^ 0 
and any n > 2. Circles of radius r = R, i.e. equators, such as the one above, are spherical choreographies for odd 
n only; for even n, it would lead to antipodal singularities. The dots show the bodies at time t — 0. 


Again, the residual of equation can be computed in Chebfun with chebop. 

As we mentioned in the introduction, the only non-circular spherical choreographies with unit 
masses found so far are for the two-body problem. Diacu and its collaborators |9j, and Perez- 
Chavela and Reyes-Victoria dzj also characterized the solutions of the spherical n-body problem, 
n > 2, in which the bodies move along the same circle (such as Figure l5T]h or along different 
ones—the relative equilibria. 

We present now new non-circular spherical choreographies. The first one is the spherical figure- 
eight , a solution of the three-body problem on the sphere of radius R = 1.4, shown in Figure [5721 
Table (5T] shows some numbers pertaining to its computation. Combining the BFGS algorithm with 
Newton’s method leads to thirteen digits of accuracy. We plot the geometrically decaying Fourier 
coefficients of the outputs of BFGS and Newton’s method in Figure [5T3l After BFGS, they decay to 
10 -7 , and after two iterations of Newton’s method, they decay to machine precision. Numerically, 
we found that the (27r-periodic) spherical figure-eight exists on spheres of radius R > 1.32. Below 
this value, it cannot fit in a single hemisphere and would therefore lead to antipodal singularities H. 

Many new spherical choreographies can be found with our algorithm. We show in Figure EH 
three spherical choreographies of the five-body problem on the sphere of radius 2. These are curved 
versions of the choreographies of Figure 13.41 Table 15.21 shows some numbers pertaining to their 
computations. We get two to five digits of accuracy with the BFGS algorithm, and applying 
Newton’s method with the outputs of BFGS as initial guesses leads to nine to thirteen digits of 
accuracy. 

Relative spherical choreographies can also be computed with this method. We plot three relative 
spherical choreographies of the seven-body problem on the sphere of radius 2.5 in Figure [53] These 
are curved versions of the relative choreographies of Figure [331 

An interactive tool to compute spherical choreographies using hand-drawn initial guesses, 


3 The solution of Figure [5~2l exists for radii R < 1.32 but with a shorter period T, with R 3 /T 2 = (1.4) 3 /(27r) 2 . 
This is a consequence of the scaling invariance described below equation 
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Fig. 5.2. Spherical figure-eight on the sphere of radius 1.4 (left) and its projection in the plane (right). The 
dots show the bodies at time t = 0. 


10 


-20 


Fourier coefficients 



Wave number 


Fig. 5.3. Absolute values of the Fourier coefficients of the spherical figure-eight of Figure f5~~2l obtained by 
BFGS (red dots) and BFGS followed by one step of Newton’s method (black circles). 



BFGS 

Newton 

Action 

18.948304138530286 

18.948304135957898 

Number of coefficients 

55 

195 

Computer time (s) 

0.51 

3.19 

Number of iterations 

72 

2 

Relative 2-norm of the gradient 

5.18e-07 

1.10e-13 

Smallest coefficient 

2.48e-07 

1.06e-17 

Relative 2-norm of the residual 

4.08e-04 

8.82e-13 


Table 5.1 

Computation of the spherical figure-eight choreography of Figure 15.21 
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Fig. 5.4. Spherical choreographies of the five-body problem on the sphere of radius 2, analogous to the planar 
choreographies of Figure \3A[ with action 58.6829 (left), 61.8035 (center), and 67.4127 (right). 



BFGS 

Newton 

BFGS 

Newton 

BFGS 

Newton 

Action 

58.6831 

58.6829 

61.8035 

61.8035 

67.4127 

67.4127 

Number of coefficients 

75 

375 

75 

205 

75 

245 

Computer time (s) 

2.35 

33.13 

1.05 

6.08 

2.14 

9.44 

Number of iterations 

142 

6 

65 

3 

115 

4 

Relative 2-norm of the gradient 

3.58e-05 

9.81e-ll 

3.50e-07 

7.61e-14 

1.98e-07 

5.72e-14 

Smallest coefficient 

1.88e-05 

2.87e-14 

4.52e-08 

1.16e-17 

6.37e-06 

2.49e-15 

Relative 2-norm of the residual 

4.26e-02 

1.68e-09 

5.62e-05 

7.32e-13 

7.05e-03 

9.32e-09 


Table 5.2 

Computation of the spherical choreographies of Figure 15.41 



Fig. 5.5. Relative spherical choreographies of the seven-body problem on the sphere of radius 2.5 with angular 
velocities 2.8 (left), —2.9 (center) and 2.31 (right), analogous to the relative choreographies of Figure 13.51 Again, 
they can be computed to ten digits of accuracy with about four hundred Fourier coefficients. 
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COMPUTING PLANAR AND SPHERICAL CHOREOGRAPHIES 



Fig. 6.1. Spherical figure-eight of Figure f5~~2l (multiplied by a factor 2) for different values of R, together with 
its planar analogue of Figure I3.lt As R increases, the spherical figure-eight converges to the planar one. 


- R = 2.0 

- R = 3.0 

R = 5.0 
— planar 





Fig. 6.2. Spherical choreographies of Figure PTD (multiplied by a factor 2) for different values of R, together 
with their planar analogues of Fiaure \3A\ As R increases, the spherical choreographies converge to the planar ones. 


choreosphere, is also available at the web-page previously given; it uses the actiongradevalsphere 
and gradhessevalsphere functions, which compute the action and the gradient, and the gradient 
and the Hessian. 

6. Limit of infinitely large radius. As its radius R gets bigger, the sphere gets flatter, and 
in the limit R oo, it converges to the complex plane. Equivalently, the spherical plane converges 
to the complex plane. The distances (|5.3[) and (|5.4[) converge to twice the absolute value, and the 
action on the sphere (|5.6[) converges to four times the action in the plane (|2.8|) . since it involves 
squares of distances. We might then expect that twice the spherical choreographies converge to the 
planar choreographies as R -A oo, and it is indeed the caseQ In Figures l6Tl 16.21 and 16.31 we plot 
the spherical choreographies of Figures 15.21 15.41 and 15.51 (multiplied by a factor 2) for increasing 
values of R and plot them together with their planar analogues. Tables l6Tl 16.21 and 16.31 report the 
oo-norm of the difference between analogous spherical and planar choreographies as R increases. It 
is clear from the tables that spherical choreographies converge to their planar analogues at a rate 
proportional to the curvature 1/R 2 . 

7. Conclusions. Choreographies are very special solutions of the n-body problem. They are 
not only periodic but also share a single orbit. We have shown in this paper that choreographies 
exist on a sphere in a cotangent potential for various n > 2. Curved versions of Simo’s planar 
choreographies, they can be computed to high accuracy using stereographic projection, trigonomet¬ 
ric interpolation, and minimization of the action. 


4 We are studying the convergence with a fixed period 27r. Similarly, T-periodic spherical choreographies converge 
to T-periodic planar choreographies. 
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R = 2.5 

- R = 3.0 

R = 5.0 
-planar 




-2 0 2 



Fig. 6.3. Spherical choreographies of Figure 15.51 (multivlied by a factor 2) for different values of R, together 
with their planar analogues of Figure 13.51 As R increases, the spherical choreographies converge to the planar ones. 



R = 1.4 

3 

5 

10 

100 

1000 


4.11e-01 

4.87e-02 

1.65e-02 

4.03e-03 

4.00e-05 

4.03e-07 


Table 6.1 

Convergence of the spherical figure-eight of Figure to the planar figure-eight of Figure EH as R increases. 


Stability properties of the spherical choreographies have not been discussed. In the plane, the 
only non-circular stable choreography is the figure-eight of Figure 13.11 We have found numerical 
evidence that the spherical figure-eight of Figure [52] is stable too. We have solved the curved 3-body 
problem (14. ip , with initial conditions defined by the reds dots (positions) and the tangents at these 
dots (velocities) of Figure [521 We ran it for a thousand full orbits, i.e., from t = 0 to t = 20007T, 
and the solution did not fall apart. All the other spherical choreographies presented in this paper 
fell apart after only a few full orbits. The systematic approach to study the stability of periodic 
solutions of dynamical systems is to compute the eigenvalues of the derivatives of the associated 
Poincare maps. We are currently working on a different algorithm, based on the singular value 
decomposition of the operator which governs the first variational equation of (|4.1D , to compute 
these eigenvalues. Details will be reported elsewhere. 

Acknowledgements. We thank Coralia Cartis and Jared Aurentz for helpful suggestions 
about quasi-Newton and Newton methods, and Alain Chenciner and Carles Simo for giving us 
details about the computation of planar choreographies. We are grateful to the reviewers for their 
comments. The first author is much indebted to supervisor Nick Trefethen for his continual support 
and encouragement. 

Appendix A. Closed-form expressions for the gradient and the Hessian in the plane. 

Let N be an odd number, and let pw(t) be the trigonometric interpolant of q{t) at N equispaced 
points on [0, 2tt) defined by ([3.1I) - ([3.2D . We can decompose the action (|2.8j) into the sum of two 
terms Ak and Ajj with q{t) and q'(t) approximated by pw(^) and p' N (t ), 

✓>277 ^ 1 /*! 

A k = — / I p'n^) +iwpN{t) I dt , Au = — I 

z J o z j =1 J o 

The two terms Ak and Au depend on the 2 N variables {uk, Vk}, \k\ < (N — 1)/2, where c& = Uk+ivk 
are the Fourier coefficients cm of PN(t )• Let V denote the gradient with respect to the u\f s and 


p N {t) Pn \ t T — 


dt. 


(7.1) 
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R= 2 

3 

5 

10 

100 

1000 

Left 

3.27e-01 

1.05e-01 

3.37e-02 

8.09e-03 

7.98e-05 

8.00e-07 

Middle 

3.06e-01 

1.04e-01 

3.39e-02 

8.17e-03 

8.07e-05 

8.04e-07 

Right 

3.34e-01 

1.12e-01 

3.64e-02 

8.75e-03 

8.65e-05 

8.65e-07 




Table 6.2 




Convergence of the spherical choreographies of 

Figure 15.41 to the planar 

ones of Fiaure\3A\ as R increases. 


i? = 2.5 

3 

5 

10 

100 

1000 

Left 

2.78e-01 

1.74e-01 

5.54e-02 

1.32e-02 

1.30e-04 

1.30e-06 

Middle 

2.19e-01 

1.39e-01 

4.50e-02 

1.08e-02 

1.07e-04 

1.12e-06 

Right 

1.69e+00 

9.47e-01 

2.35e-01 

5.22e-02 

5.04e-04 

5.02e-06 


Table 6.3 

Convergence of the spherical choreographies of Figure 1531 to the planar ones of Figure 1531 as R increases. 


f/c’s, that is V = (V u ,V v ) T , V u = (d/du k ) T , V v = (d/dv k ) T , |fc| < (N — l)/2. We wish to derive 
closed-form expressions for W Ax and V 'Ajj. Consider first with 


N-l 

2 


N-l 

2 


A K(u k ,v k )=Trn ^2 \(k + u>)c k \ 2 = im ^ (k + w) 2 (u| + v\), 


(7.2) 


k=-- 


k=— - 


since p' N (t ) + iujpx(t) has Fourier coefficients {ikc k + iujcj . c }, and using Parseval’s identity. It leads 
to 


du k 

Consider now Ajj, with 

n-l «2 

Au{u k ,v k ) = — V] / 

2 .7 — 1 


3 Ax /7 xo 3 Ax . . o . . TV — 1 

= 2nn(k + uj) 2 u k , ——= 2nn(k + uj) 2 v k , \k\ <-. 


<%/c 


dt 


= ? y fj(u k ,v k ,t) = Pat(0 


J=1 - ~ V fj(u k ,v k ,t) 

Expanding PN(t) and + 2irj/n) and regrouping real and imaginary parts lead to 

/ ^ \ 2 / ^ \ 2 
fj(u k ,v k ,t ) = a kj(t) u k + b kij (t)v k + 5^ a k j(t)vk ~ b k j(t)u k 


<k=- 


< k=-~ 


with 


ctkj(t ) = [l — cos(27rjA:/n)] cos(fct) + sin(27r//c/n) sin(fct), 
b k j(t) = [ — 1 + cos(27r/&/n)] sin(fct) + sin(27r/£;/n) cos(fct), 


(7.3) 


(7.4) 


(7.5) 


(7.6) 


for |fc| < (AT — l)/2 and 1 < j < n — 1. The partial derivatives of Ajj with respect to the u k s and 
v k s can then be computed with the chain rule, 


n ^ r 2n 

VAu =2^i V 

j=l 


l 


vT 


n-l ,277 


dt — — 


II — X 


V/i 


=1 / 


3/2 


dt , 


(7.7) 
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with 


df: 


7 ^ = 2 ^2 ( a i,j(i) a k,j(t) + bij(t)b k ,j(t) u t + bij(t)a kj (t) - aij(t)b k j(t) vi ), (7.8) 


i=-i 


and 


df± 

dvk 


N-1 
2 


= 2 a u( t ) b k,j( t ) ~ b ij(t) a kj(t) ui + bij(t)b kj (t) + aij(t)a kj (t) vi . (7.9) 


i=-± 


Let us now derive the formula for the exact Hessian matrix H. 

' dAA _jPA_ 1 

duidu k duidv k 

H = 

' (PA d\A 

dviduk dvidvk - 


(7.10) 


H is a 2 N x 2 N matrix, and each block is N x N. Note that the d 2 A/dviduk block is the transpose 
of the d 2 A/duidvk block, so we are going to derive formulas for the d 2 A/duiduk, d 2 A/duidv kl and 
d 2 A/dvidvk derivatives only. It is clear from (|7.3|) that 


d 2 A k d 2 A k 2 d 2 A k 

^— = o ^ = 27 rnk dki, ~ ~ 

duiduk dvidvk ouidv k 


< 


N -1 


(7.11) 


where S k i is the Kronecker delta. 

The second derivatives of Ajj with respect to the u k s and v k s can be obtained by differentiating 
G2D one more time, e.g., 


iPAu 

duidu k 


-i Si 


-1 /* 2 tt 9 fj f _ 3 dfj dfj 

duiduk ■'•7 2 dui duk 


, 5/2 


dt , 


with 


d 2 fj 


— 2 (aijdkj + bijbkj). 


duiduk 

There are similar formulas for the other derivatives with 

0 2 f) _ d*p d 2 fj 


dvidvk duiduk duidv k 


— 2 ( a i,jbk,j bija k ,j). 


(7.12) 


(7.13) 


(7.14) 


Note that all the second derivatives involving the real and imaginary parts uo and vo of the 
constant terms cq are zeros, i.e., 


8\Ak _ <PAk _ <PAk 
duiduo dvidvo duidvo 


(7.15) 
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and 

d 2 Ajj d 2 Ajj _ d 2 Ajj _ g 

duiduo dvidv o duidv o 

To prove (|7.15p . take k = 0 in (JTTTTj) , and to prove (|7.16p . note that clqj = boj = 0. As a 
consequence, when using Newton’s method with the exact Hessian (Pol) . one needs to get rid of 
these derivatives to make the matrix nonsingular, that is, eliminate the lines and columns that 
correspond to the constant term. Similarly, one needs to get rid of uq and vq in the vector of 
optimization variables. 

Appendix B. Closed-form expressions for the gradient and the Hessian on the 
sphere. Again, let N be an odd number, and let pN(t) be the trigonometric interpolant of q{t) at 
N equispaced points on [0, 2tt) defined by (|3.ip ~ (j3.2p . We can decompose the action (|5.6p into the 
sum of two terms Ak and Ajj. The first term comes from the kinetic energy, 


A K 


n f 27T / 2R 2 \p' N (t) + iup N (t )\\ 2 

2 J 0 \ R 2 + \PN(t)\ 2 J 


while the second term comes from the potential energy, 



(7.17) 


(7.18) 


with 




2R 2 -DAt) 2 

- ; J = l =, 1 < ? < n — 1 


(7.19) 


with Dj(t) = d(p N (t),p N (t + )). 

Again, let V denote the gradient with respect to the s and f/ds, that is V = (V w ,V v ) T , 
V w = {d/duk) T , V v = (d/dvk) T i \k\ < (N — l)/2. Let us first derive the closed-form expression 
for WAk- A straightforward calculation leads to 




ru 


2-R 2 |Pjy (t) + uop N {t) [ \ 

R 2 + \p N {t)\ 2 J 


2-i 


dt = 2 nR‘ hVs :, sVh dt, (7.20) 

Jo 


h 2 


with 


g(u k ,v k ,t) = | p' N (t) + iujp N (t) | 2 , h(u k ,v k ,t ) = (i? 2 + \pN(t)\ 2 ) . 
The functions g and h are given by 

JV — 1 2 i— N — 1 


9 = 


and 


(k+u)(uksin.(kt)+VkCos 


L k=— - 


+ 


(fc+o;) (uk cos (kt) —Vk sin (kt)^ 


k —— - 


(7.21) 


, (7.22) 


h = 


A> 2 + £ Uk sin(H) + Vk cos(kt) j + £ t* cos(Jfc) - sin(Atf) 




^ k=- + 


(7.23) 
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Their partial derivatives are given by the formulas 

N-l 


= 2 (k + id) ^ (Z + cj) fui cos ((fc — Z)t) + vi sin ((fc — 
O'U'k -- - ' 


z=-^ 


-^- = 2 (k + u>) ^2 (l + u)(ui sin ((l - k)t) + Vicos ((l - k)t)^), 


i=-i 


and 


5h 


— =aVh E ui cos ((fc — Z)t) + vi sin ((fc — Z)t), 




z=-^ 


= 4V7z ui sin ((/ — fc)t) + vi cos ((/ — fc)t). 


i=- + 


Let us now derive the closed-form expression for VAjj, 

n-l r 2ir 


f‘ Z7T 

VA u = E l vAjjityit. 


with 


Let us write 


v4w = 


-8 R 4 VDj(t) 


UyU) - „ ^N 2 l3/2 


^i(i) = 


Dj(t)*[4B? - Dj{t)* 
2 J R 2 |p7v(t) — PN{t + ^p)| 


, 1 < j < n VI. 


= 2i7 


2 V77 


(7.24) 


(7.25) 


'(i?2 + | Piv (i)| 2 ) (i?2 + | pAr(t + 23i) | 2 ) ^ 

with fj{uk, Vk , £) = IpivW — Pa(^ + 27rj/n)| as defined in (|7.5D - (|7.6|) . and defined by 

rj(u k , v k ,t) = \JrA- p N (t + ^) , 


(7.26) 


(7.27) 


(7.28) 


(7.29) 


that 


Ti = 


\ 


T I ^ ^ Ck,j (fy'U'k T dk,j(t)vk I T | ^ ^ Ck,j dkj (i'jujz J , (7.30) 


v k=~- 


< k=— - 


with 


Ckj(t ) = cos(27t kj/n) cos (kt) — sin(27rhjf /n) sin(fct), 
dk,j(t) = — cos(27t kj/n) sin (kt) — sin(27rhj/n) cos (kt), 


(7.31) 
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for \k\ < (N — l)/2 and 0 < j < n — 1. It leads to 


VD At)= 2R ^ r ’ Vf ‘ / 2 2 -fL V(V ‘> , 1 

r o r j V Jj 


< j < n — 1. 


(7.32) 


The derivatives of fj with respect to the u k s and v k s are given by (|7.8[) - (|7.9p . while the derivatives 
of rj are given by 


drj 1 
8u k Tj 


c k,j (t) Cl ’i W Ul + d k3 (f) V l - d kJ (t) ^ vi ~ d k3 (t) U l 


<l=- + 


>l=- + 


dv k rj 


TV — 1 
2 


TV-1 

2 


— - — { y] Cij{t)ui +di t j(t)vi ) +c k j(t)[ y] Cij{t)vi - dij(t)ui 


'/=-■! 


W=--i 


(7.33) 


for |fc| < (N — l)/2 and 0 < j < n — 1. 

Let us now derive the formula for the exact Hessian H on the sphere. Differentiating (|7.20[) 
gives the second derivatives of Ak with respect to the u k s and v k s, e.g., 


d 2 A 


K 


duidu k 


_ 2 |TLp 

= 2ni? 4 / ^- ,* 9ui dt, P k = h 


'L 


h 3 


dg dh 
du k 9 du k ' 


with 


and 


<% dg + ^ dh d 2 h 


d 2 g 


dui dui du k duiduk dui du k ^ duidu k ’ 

d 2 h 


= 2(fc + cj)(Z + lS) cos ((fc — Z)t), 


duiduk v duiduk 2h dui 

There are similar formulas for the other derivatives, with 
d 2 g d 2 g d 2 g 


= 4lr+ 4v ^ cos (( fc - oo- 


dvi dv k dui du k ’ dv k 


= 2(fc + u;)(Z + cj) sin ((/ — fc)t), 


and 


d 2 h 1 dh r- f 7X v d 2 h 1 dh r- . , 7N v 

+ 4vW ((fc - 00, 777777^7 = ^ 77777 + 4v ^ sm U “ *)*)• 


dv\dv k 2 h dvi v } duidv k 2 h dui 

The second derivatives of can be obtained by differentiating (|7.27D , e.g., 




16T? 4 ^f^- 

an; an fc 


241? 4 

an; oufe 


8i?‘ 


l4 _aAD 2 _ 

duidu k 


duidu k d?[ 4R 2 -Dff /2 Dj [4R 2 - D 2 ] 5/2 £>? [4i? 2 - £)?] 3/2 ’ 


with 


5 2 D 


duidu k 


j _ 2 R 2 rorjj At Z d_A irof iM + 2r Uiii7 + Po r i§ft) 


r 3 r 3 f ^/2 

r 0 Vi 


(7.34) 


(7.35) 


(7.36) 


(7.37) 


(7.38) 


(7.39) 


(7.40) 
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and 


1 dfj t dr 0 drj 

qk = 2 v ‘a^- , ' r '^ !in d^' 


(7.41) 


The second derivatives *rm involve the derivatives of can with respect of the ui s-—the chain 
rule leads to nine terms, and only two of them are new, d 2 r^jdu\duk and d 2 rj/duiduk- These are 
given by 


duidu k 


c k,j c i,j + dkjdij Qy^ 


0 < j < n - 1. 


(7.42) 


Similar formulas can be derived for the other derivatives, with 


d 2 rj _ c k,jQ,j + dh,jdij dv 3 k d 2 rj _ c i,jdk,j Ckjdij du 3 k 

dvidvk rj ’ duidvk rj 


(7.43) 


Note that, as in the plane, all the second derivatives involving the real and imaginary parts uq 
and vq of the constant terms Co are zeros, so when using Newton’s method with the exact Hessian, 
one needs to get rid of these derivatives to make the matrix nonsingular. 
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